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Abstract. A pseudo-Newtonian description of the grav- 
itational field which yields the equations of motion re- 
sembling, as closely as possible, the geodesic equation of 
general relativity is found for the Kerr spacetime. The 
potential obtained consists of three parts, interpreted as a 
purely Newtonian (gravitoelectric) term, a dragging (grav- 
itomagnetic) term, and a space-geometry correction. The 
accuracy of the pseudo-Newtonian model is studied by a 
method which compares, systematically, two large sets of 
trajectories: geodesies in the Kerr spacetime versus test- 
particle trajectories in the pseudo-Newtonian field. It is 
suggested that every pseudo-Newtonian model should be 
submitted to analogous systematic analysis before it is 
used in astrophysical applications. A modified Newtonian 
potential which accounts for the frame-dragging effects 
can be a practical tool in studying stationary accretion 
discs. Non-stationary configurations are more complicated 
and we will not suggest to use this approach to such topics 
in accretion theory as gravitomagnetic oscillations of discs 
and their relation to quasi-periodic sources. 

Key words: Accretion: accretion-discs - black hole 
physics 



1. Introduction 

Most of theoretical astronomy uses Newtonian theory of 
gravitation, considering that the effects of general relativ- 
ity are weak for most astronomical objects. Even in sit- 
uations where these effects do become important or even 
dominant, pseudo-Newtonian models have often been ap- 
plied which aspire to mimic the corresponding relativistic 
situations. Pseudo-Newtonian approaches aim to save nu- 
merical work, to enable analytical solutions where fully 
relativistic treatment is too cumbersome, and to provide 
new insights into the relativity theory and its implications. 



The present contribution concerns pseudo-Newtonian 
models developed in order to describe very compact ob- 
jects which presumably reside in the nuclei of galaxies 
and in some X-ray binaries (e.g., Rees 1998). In these 
sources, a key feature is an accretion disc around a ro- 
tating black hole. The actual accretion flows are most 
likely non-stationary, non-axisymmetric, and described by 
complex local physics. For the sake of simplicity, how- 
ever, a standard model of disc accretion (see Kato, et 
al. (1998) for a textbook exposition of the accretion the- 
ory including its recent advances) provides an ingenious 
approximation: the standard disc is smooth, axially sym- 
metric, geometrically thin, and characterized by three pa- 
rameters. It is astrophysically realistic only in a restricted 
range of its parameters. The pseudo-Newtonian approach 
has been devised in order to introduce effects of gen- 
eral relativity into accretion models, especially in the case 
of geometrically thick discs (Paczyhski & Wiita 1980). 
The black hole is treated there as a Newtonian body 
whose gravitational field is determined by the potential 
V pw = —M/(r — 2M) (geometrized units with c = G = 1 
will be used throughout the paper). This simple expression 
mimics the gravitational field of a non-rotating (Schw- 
arzschild) black hole quite accurately, in particular, it re- 
produces correctly the radii of the marginally stable and 
the marginally bound circular orbits of free test particles, 
r ms = 6M and r mb = AM, Recently, Abramowicz et al. 
(1996) proposed a certain rescaling of velocities in the field 
of V pw to obtain better agreement with corresponding rel- 
ativistic values. In this way, calculations of the observed 
spectra of the discs could also be improved. It has been 
argued, however, both on theoretical grounds (Bardeen 
1970) and from observations (Iwasawa et al. 1996; Karas 
& Kraus 1996), that the central black holes in galactic nu- 
clei would be more likely to rotate rapidly. Thus the main 
feature that should be included in the pseudo-Newtonian 
models is the rotation of the central object. In addition, 
the rotation-induced dragging must be taken into account 
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when studying non-equatorial warped discs around com- 
pact objects (Bardeen & Petterson 1975). 

The pseudo-Newtonian potential for a non-rotating 
black hole has been used in numerous works and we recall 
at least a few of them. For example, Abramowicz et al. 
(1988) used V pw in their introductory work on slim discs. 
Okazaki et al. (1987) adopted this approach to describe 
global trapped oscillations of relativistic accretion discs. 
Later, Nowak & Wagoner (1991) devised another form of 
the potential which is better suited for their purpose be- 
cause it reproduces the epicyclic frequency k with higher 
accuracy than V pw . Szuszkiewicz & Miller (1998) studied 
the limit-cycle behaviour of thermally unstable flows, also 
with the help of the pseudo-Newtonian description. On 
quite a different subject, Daigne & Mochkovitch (1997) 
applied V pw to study the runaway instability in accre- 
tion discs which could trigger gamma ray bursts. Very 
recently, Ruffcrt & Janka (1998) used V pw in a detailed 
study of neutron-star mergers and related production of 
the bursts. It is quite apparent from this short list of dif- 
ferent astrophysical applications that the approach has 
its advantages and enables qualitative (and quantitative, 
with some caution) studies of different topics concerning 
black holes. Especially the problems of accretion can be 
treated in this manner. It is understood that quantitative 
results need always be checked carefully within the full 
relativistic framework. This is particularly true for time- 
dependent phenomena (like oscillations and waves in flu- 
ids) and for computations of observed spectra. Notice that 
a textbook overview of the whole subject can be found in 
Kato et al. (1998). 

Now we come to the case of a rotating black hole. At 
least two topics of current astrophysical interest can im- 
mediately be mentioned where the frame-dragging effects 
around a rotating black hole are important: the problem 
of oscillations of relativistic accretion discs (see Wagoner 
1998 for a recent review), and that of precessing discs in 
low-mass X-ray binaries (Stella & Vietri 1997) and black- 
hole binaries (Wei et al. 1998). It is apparent that the 
whole subject of oscillation modes in relativistic discs calls 
for detailed investigation (cf. also Markovic & Lamb 1998) 
and a suitable pseudo-Newtonian formulation can be an 
appropriate tool before embarking on a complete solution. 
However, unlike the case of stationary gaseous configu- 
rations, nonstationary phenomena require a more refined 
choice of the pseudo-Newtonian model because there are 
additional quantities apart from the location of marginal 
circular orbits which must be correctly modelled (e.g. n). 

Although the ever increasing computational facilities 
will make practical reasons for the pseudo-Newtonian ap- 
proach rather old-fashioned in near future, what still re- 
mains desirable is a simple expression which would sim- 
ulate the rotating (Kerr) black hole accurately. This is 
however difficult to find. For example, the potential sug- 
gested by Chakrabarti & Khanna (1992) could be useful in 
studying thin accretion discs around rotating holes, but it 



applies only to the equatorial plane and its interpretation 
is rather unclear (free parameters are fitted in a purely 
pragmatic manner for a restricted set of trajectories). An- 
other form of the pseudo-Newtonian potential has recently 
been proposed by Artemova et al. (1996). This pseudo- 
potential reproduces very well steady circular motion but 
it ignores all effects which in the true Kerr metric are as- 
cribed to the Lense-Thirring precession and which make 
test-particle trajectories non-planar. This is also the main 
reason why in contrast to the non-rotating case, pseudo- 
potentials for the Kerr metric have had rather restricted 
impact. Also, a proper understanding of these generaliza- 
tions may be as difficult as employing general relativity 
fully. However, the motivation which stems from attempts 
to understand and interpret predictions of the relativity 
theory has not disappeared. 

It is the aim of our present contribution to discuss the 
pseudo-Newtonian modelling of the gravitational field of a 
rotating (Kerr) black hole, and also to propose how to test 
such models. One should remember that the idea of the 
pseudo-Newtonian approach represents a certain mathe- 
matical model rather than a rigorously defined approxi- 
mation such as e.g. the weak-field approximation of the 
Einstein equations. The model does not aspire to repre- 
sent any gravitational theory and, in particular, the cor- 
responding potential is not required to satisfy any field 
equations. (It is exactly the lack of precise definition of 
the approximation method which leads us to introduce 
a test of accuracy of the model in the present article.) 
This fact, however, does not diminish the practical value 
of Paczynski-Wiita potential and other pseudo-Newtonian 
models, which apply even to regions with strong gravity 
and capture qualitative features of motion near the hori- 
zon. 

We start by writing down the spatial components of 
the Kerr geodesic equation in Boyer-Lindquist coordinates 

= (t,r,0,# 

Sr = MAXT 2 (£ - 2r 2 )(t- a</>sin 2 <9) 2 

+[(r - M)S/A - r] f 2 + a 2 r9 sin 26 

+rA(6 2 + 4> 2 sin 2 9) , (1) 

W = {2Mr^- 2 [ai - (r 2 + a 2 )4>] 2 + a 2 2 - 

a 2 r 2 /A + A0 2 } sin (9 cos (9- 2rf9, (2) 
±A£ 2 = Mo(E -2r 2 )(i- acj> sin 2 6)r 

+2MrA[ai - (r 2 + a 2 )</>] 9 cot 9 

-£(£- 2Mr)(rr + A<9 cot 6>)</> , (3) 

where 

A = r 2 -2Mr + a 2 , S = r 2 + a 2 cos 2 9; (4) 

M and a are parameters of the Kerr solution, and the dot 
denotes differentiation with respect to the affine parame- 
ter normalized so that the 4- momentum is = x^. 

In the next section, Newtonian equations of test- 
particle motion in an axially symmetric gravitational field 
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are given. In Sect. we derive a pseudo-Newtonian poten- 
tial for which these equations get a form very similar to the 
above geodesic equation in the Kerr spacetime.EJ In Sect. 
[|, numerical integration is carried out for a large number 
of trajectories both in the Kerr and in the "pseudo-Kerr" 
fields. A specific criterion analogous to the method of Lya- 
punov coefficients is introduced and the rate of divergence 
of each couple of corresponding trajectories is determined 
numerically. The mean values of this rate obtained for 
various combinations of constants of motion indicate the 
quality (i.e. accuracy, as defined below) of the pseudo- 
Newtonian approximation. 

2. Stationary and axisymmetric field 

In Euclidean space and in ellipsoidal coordinates x l = 
(r, 0, 0), related to the Cartesian coordinates (x, y, z) by 



3. Potential for the Kerr field 

Relativistic gravitational fields can be interpreted as con- 
sisting of three parts: the Newtonian (also Coulomb or 
gravitoelectric) component, the dragging (gravitomag- 
netic) component, and the space-geometry component. 
The Newtonian component is generated by mass-density 
and is given essentially by the gradient of the g 00 com- 
ponent of the spacetime metric. The dragging compo- 
nent is generated by mass-currents and determined by 
curl of <7oi- The space-geometry component (determined 
by gij) has no classical analog. This view and, in particu- 
lar, the analogy with electromagnetism are most straight- 
forward within the linearized approximation, but can be 
given a precise and invariant meaning even in strong fields 
(Thorne et al. 1986; Jantzen et al. 1992, and references 
cited therein). Let us compose the potential for the Kerr 
spacetime in a way that acknowledges this approach. 



x = v r 2 + a 2 sin0cos</>, 
y = \J r 2 + a 2 sin sin < 



= r cos 8, 



the kinetic energy of a particle of mass m is 

«2 



T 



P 
2m 
1 
2m 



x 2 + y 2 + z 2 



2m 

^2 



+ 2 +(r z +a')^sm z 



(5) 



where a is a non-negative constant of the dimension of 
length and the dot denotes differentiation with respect to 
the time parameter t/m (normalized so that the momen- 
tum is p = x). The motion of the particle in a stationary 
axisymmetric potential U = U(r, 8, r, 8, <j>) is described by 
equations 



Sr 



m 2 (r 2 + a 2 )[(U t r)' - U. r ] 

-r(af sin 8) 2 /{r 2 + a 2 ) + a 2 r8 sin 28 

+r (r 2 + a 2 )(0 2 + </> 2 sin 2 8), (6) 



£0 = m 2 [(U ( 



U.t 



2 -2 ;/ 2 i 2\ 

a r /{r + a ) 



2 ] sin cos - 2rfd , 



(r 2 + a 2 ) 



m(U^Y sin -2 8 

-2[rr+ (r 2 +a 2 )0cot 8](f>, 



where we denote, for example, 

^ =m dt u 



(7) 
(8) 

(9) 



For a static and spherically symmetric field (both gravita- 
tional and non-gravitational) , a similar problem was solved by 
Jaen (1967) using the Hamilton- Jacobi theory. 



3.1. The Newtonian component 

Cutting the Kerr manifold at r = 0, one can obtain the 
spacetime which is free of causality- violating regions. The 
cut leads to jumps in derivatives of the metric at r = 
which are interpreted as a thin massive layer. This induced 
mass spreads over the r = hypersurface and consists of 
the attractive singular ring (z = 0, p f = \/ x 2 + y 2 = a) 
of infinite positive mass density, spanned by the repulsive 
disc (z — 0, < p t < a) of negative mass density. Con- 
structing the causally maximal extension of the Kerr met- 
ric, Keres (1967) and Israel (1970) showed that the New- 
tonian field generated by the massive layer corresponds to 
the potentiala 



V = -Mr/S. 



(10) 



The structure of this field was shown in Semerak (1995; 
cf. Fig. 2 there) by depicting its lines in the (p t , z)-plane. 
We will use the Keres-Israel potential as a scalar "seed" 
of our model. With U = V given by Eq. @, Eqs. (6)-(8) 
read 



m 2 M (r 2 + a 2 )(S - 2r 2 )/E 2 

-r(af sin0) 2 /(r 2 + a 2 ) + a 2 r8 sin 20 



+r(r 



a 2 )(0 2 



(j) 2 sin 2 i 



£0 = {2m 2 Ma 2 r/Y, 2 + a 2 8 2 - a 2 f 2 /{r 2 + a 2 ) 
+ (r 2 +a 2 )</> 2 }sin0cos0- 2rf0 , 
$ = -2 [rr/(r 2 + a 2 ) +8 cot 8]4>. 



(11) 

(12) 
(13) 



Separated first integrals of these equations (analogy of 
Carter's equations for the Kerr spacetime) were found by 



2 Cf. also Qadir (1986) where an electrically charged gen- 
eralization of this potential was discussed. A somewhat dif- 
ferent expression, V = —(M/a) arctan(a/r), was obtained by 
Krasiriski (1980) in construction of a Newtonian model of the 
Kerr source. For the other result, the potential of an oblate 
spheroidal homeoid V = — (M/a)arccot(r/a), see Misra (1970). 



4 



O. Semerak & V. Karas: Pseudo-Newtonian models of black holes 



Israel (1970). It was illustrated in Semerak (1996) that 
Eqs. (11)-(13) often yield trajectories indiscernible from 
their exact-Kerr counterparts computed from Eqs. (1)- 
(3). However, they do not contain the terms linear in ve- 
locities which embody the very characteristic feature of 
the Kerr geometry — the frame-dragging effects. 

3.2. The dragging component 

In this section, we are led by analogy with classical elec- 
trodynamics and by observation that the Kerr dipole-like 
gravitomagnetic field resembles the Kerr-Newman mag- 
netic field. One thus arrives at the potential which incor- 
porates dragging in the form 



U = V-v-A, 



(14) 



where v = p/m and the vector potential A is given by 
that of the Kerr-Newman electromagnetic field, Aj = 
(0, 0, Qrasin 2 #/£), with Q (electric charge of the Kerr- 
Newman centre) replaced by — 2M (the extra factor of 
2 is ascribed to the tensorial character of gravity and 
the minus sign to its attractive nature) — i.e., Ai = 
(O,O,2F<zsin 2 0) = (0,0, g^). Similar (just half) correc- 
tion for the dragging was considered by Dadhich (1985) 
for a special case of particles with zero axial angular mo- 
mentum. The added term — v ■ A really introduces the 
desired dragging terms into the Eqs. (11)— (13): they read 
now 

Sf = [r.h.s. of Eq. (11)] - 2mMT,~ 2 

x(r 2 + a 2 )(£-2r 2 )a0sin 2 6>, (15) 
£0 = [r.h.s. of Eq. (12)] 

-2mMrYT 2 (r 2 +a 2 )a0 sin 20, (16) 
i(r 2 +a 2 )S 2 ^ = mMa[(Y,-2r 2 )r 

+2r(r 2 + a 2 ) (9 cot 0] 

-£ 2 [rr + (r 2 + a 2 ) cot 0]4>. (17) 

Note that the contravariant 3-potential obtained by using 
the inversion of the Kerr 3-metric, g %J — \/g%j, is 



A i = (0, 0, ga/gto) = (0, 0, -w K ) > 



(18) 



where uj k = 2Mar/A and A = AS + 2Mr(r 2 + a 2 ); A 1 
equals the shift vector of Thorne et al. (1986) which stands 
for the potential of the gravitomagnetic field in this refer- 
ence. 

3.3. The space- geometry component 

The space curvature cannot be understood directly within 
the Newtonian or electromagnetic analogy. It may only be 
included by introducing ad hoc corrections into the form 
(HI: 



U = V - v ■ A 



(v ■ A) 2 

AV 



MrY, 




1 a sin 

m 



(r 2 + a 2 ) A m 2 



(19) 



This potential leads to equations of motion 

Sf = MAIT 2 (£-2r 2 )(m-a0sm 2 0) 2 

+[(r- M)T,/A-r] r 2 

+a 2 r9 sin 29 + rA{9 2 + 2 sin 2 9) , 
£0 = {2MrY>- 2 [am- 

(r 2 + a 2 )0] 2 + a 2 9 2 

-a 2 r 2 /A + A<£ 2 }siri0cos0- 2rr9 , 
±ZA4> = Ma(S -2r 2 )(m-a4) sin 2 9)r 

+2Mr(r 2 + a 2 ) [am - (r 2 + a 2 )<j>] 9 cot t 

-S 2 (rf + A9 cot ( 



(20) 



(21) 



(22) 



The only point in which Eqs. (20) and (21) differ from 
the relativistic Eqs. (l)-(2) is that the relativistic variable 
i — given by (AY.)~ 1 (AE - 2Mar<&) (E and $ stand for 
the particle's energy and axial angular momentum at in- 
finity) — is replaced by m in the classical model. This 
distinction reflects, however, the conflict between the very 
roots of classical physics (where time parameter is univer- 
sal) and relativity (where proper time and some coordi- 
nate time occur, related to each other in a specific way at 
each point). Notice that Eq. (22) differs from (3) also in 
several other points (namely M is neglected three times). 

3.4- Specific features of motion in the pseudo-Kerr poten- 
tial 

Independence of the Lagrangian L = T — mil on t and 
implies two constants of motion, 



E 



mv 2 /2 — mMr/S 



1 

2m 



—d> 2 sin 2 ( 



mMr 



mA , . . 2 

— (w-w K )sw i 



(23) 
(24) 



where u> = d(f>/dt — (p/m. The Kerr axial angular mo- 
mentum at infinity contains (again, as in Eqs. (20)-(21)) 
t instead of m (u> K is the angular velocity with which the 
field rotates relative to an observer standing at infinity). 
For large r one obtains the usual forms of the energy and 
axial angular momentum in the Newtonian central gravi- 
tational field: 

E ~ (2m)" 1 (f 2 + r 2 (9 2 + r 2 2 sin 2 9) - mM/r, (25) 



$ ~ mr 2 <b sin 2 I 



(26) 



(r 2 + a 2 ) A m 2 



Let us determine the acceleration a 1 of an observer or- 
biting uniformly (oj = const) at r = const, 9 — const. 
We will define acceleration as a specific force necessary to 
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keep the observer in the orbit, i.e. as the minus accelera- 
tion of a free particle having r = 9 = at a given point. 
According to Eqs. (20)-(22), we find 

a r = AXT 3 [Af(2r 2 - S)(l - aw sin 2 6») 2 

-r(£wsin(9) 2 ] , (27) 
a e = ~S- 3 {2A/r[a - (r 2 + a 2 )uj] 2 

+AE 2 ^ 2 } sin 9 cos 9 , (28) 
a = 0. (29) 

This 3-acceleration differs from the space part of 4- 
acceleration of the Kerr stationary observer only by the 
absence of the multiplicative factor (u 1 ) 2 (square of the 
time-component of the observer's 4-velocity) [cf. Semerak 
1993, Eqs. (37)-(39)]. 

For a static observer (to = 0) one obtains, in particular, 

a' = (M/E 3 ) (A(2r 2 - S), -ra 2 sin 26, 0) . (30) 

Hence, the field is repulsive at < r < a|cos6*| in the 
sense that the radial component of j3C| ) is negative there. 

According to (27)-(29), the stationary observer needs 
no thrust (i) at r = a on the axis (9 = 0°, 180°), and (ii) in 
the equatorial plane if his angular velocity is u> = l/(o ± 
^/r 3 /M). The latter agrees exactly with the Keplerian 
angular velocity of an equatorial observer orbiting along 
a circular geodesic in the Kerr spacetime. 

Important circular geodesies in the equatorial plane, 
on the other hand, are not reproduced properly by the 
potential (19), viz. the equation for marginally stable or- 
bits reads 

r 3 ± 8aVMr 3 - 3a 2 r - 2Ma 2 = (31) 
and for marginally bound orbits 

r 2 ± 4aVMr - a 2 = 0. (32) 

The correct equations have respectively the forms 

r 2 - QMr ± SaVMr ~ 3a 2 = (33) 

and 

r 2 - 4Mr ± 4ay/M^ - a 2 = 0. (34) 

In the Schwarzschild case, for instance, both our equations 
imply r — 0. 

In any pseudo-Newtonian description of a rotating 
black hole, the main difficulty is to simulate the pres- 
ence of the horizon and of the dragging effects with an 
acceptable precision. Both these phenomena are outside 
the scope of Newtonian physics. Our potential accounts 
for the frame-dragging effects, and especially in the inter- 
mediate and large distances provides a good fit, whereas 
it does not describe correctly the innermost region where 
the horizon and important circular orbits lie. In order to 



account for the horizon, one would have to start from the 
scalar potential which diverges to — oo there. The Keres- 
Israel potential (pX|), instead, reaches — oo only at the very 
singularity E = 0. The potential proposed by Artemova et 
al. (1996) is a better alternative in this respect, being the 
simplest generalization of the Paczyhski- Wiita potential 
to the Kerr case which reproduces the horizon and approx- 
imates the important orbits. Also, the epicyclic frequency 
of small radial oscillation n(r), important in the theory 
of discoseismology, is not reproduced with good accuracy 
by potentials ([!(]) and (19) although it does show a max- 
imum, typical for a relativistic k. Apparently, trying to 
incorporate different manifestations of the frame-dragging 
accurately, one may end with a long expression without 
practical sense. 

To summarize, each particular relativistic effect can 
be well simulated within classical physics, but pseudo- 
Newtonian modelling of the complete relativistic situation 
has only a restricted validity. In order to clarify the value 
of our potential, we have carried out an extensive com- 
putation of trajectories and introduced a criterion which 
determines the accuracy of the pseudo-Newtonian model. 

4. Testing approximative equations 

As mentioned above, the pseudo-Newtonian potential has 
been used frequently in various situations in which general 
relativistic effects on the motion of material (either test 
particles or fluids) are essential but exact calculations are 
too difficult. It is however impossible to estimate, a priori, 
the error which is introduced by replacing the original sys- 
tem, described in the framework of general relativity, by 
a corresponding pseudo-Newtonian system. The plausibil- 
ity of a particular form of the simulating potential can be 
verified by solving analogous situations within the exact 
theory. For example, several specific questions in the as- 
trophysics of accretion discs had been first analysed using 
the pseudo-Newtonian theory, and the results were only 
later supported by more complicated calculations within 
the Schwarzschild spacetime (with identical local physics 
of the fluid). Indeed, it is quite trivial to recall that the 
above-mentioned standard model of thin discs was formu- 
lated within the Newtonian (Shakura & Sunyaev 1973), 
pseudo-Newtonian (Paczyhski & Wiita 1980), and also rel- 
ativistic frameworks (Novikov & Thorne 1973). A specific 
response of relativistic accretion discs to oscillations (Kato 
& Fukue 1980) was also explored with a modified Newto- 
nian potential (Nowak & Wagoner 1991) before further 
steps were carried out (Perez et al. 1997). These examples 
indicate that the pseudo-Newtonian model is appropri- 
ate for investigating the motion of material around black 
holes. 

It is quite straightforward to check whether the 
pseudo-Newtonian potential is suitable for treating the 
motion of test particles and fluids when one deals with 
a spherically symmetric system. What appears more dif- 
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ficult is to develop an analogous pseudo-Kerr theory, de- 
scribing the relativistic effects near a rotating compact 
object. Now we will propose a systematic approach to es- 
timate the quality of a particular model by computing 
a sufficiently large number of trajectories with different 
initial conditions and comparing results with exact calcu- 
lations of geodesic motion. The test we suggest determines 
the rate of divergence of trajectories (as given by approxi- 
mative versus exact equations). We propose that this type 
of check should be carried out for each particular set of 
approximative equations before they are applied to as- 
trophysical situations. (Until now, only rather restricted 
tests on a relatively small number of trajectories have been 
applied. As a consequence, one cannot be sure about pre- 
dictions based on approximative equations of such models 
because there is little or no control of their overall preci- 
sion.) For illustration, we then submit the approximative 
equations of Sect. || to our test. 

4--1- Criterion of accuracy of approximative solutions 

Our criterion is motivated by the definition of the Lya- 
punov characteristic coefficients which have been intro- 
duced in order to characterize chaotic systems (see, e.g., 
Chapt. 5 of Lichtenberg & Lieberman (1983). The most 
important features of the Lyapunov coefficients relevant 
for the present work are summarized in the Appendix. 

We will now introduce a parameter which is analo- 
gous to the maximum Lyapunov characteristic exponent 
A. One should however realize the basic difference between 
the standard formulation of the problem of chaotic system 
and our present situation. While A characterizes the rate 
of divergence of two neighbouring (initially close to each 
other) trajectories, in the present work one always starts 
with exactly identical initial conditions. Trajectories then 
separate because the motion of the first test particle is 
determined by the geodesic equation in the Kerr space- 
time, while the other (fictitious) particle moves according 
to approximative equations. The geodesic equation could 
be integrated analytically but not the approximative equa- 
tions, so we have to resort to numerical solution. Two 
points should be mentioned: (i) as the test particles move 
in a stationary system with preferred time coordinate (i), 
the separation of trajectories in the phase space is calcu- 
lated in the t — const slice; (ii) periodic rescalings of the 
separation w to zero length keep the two trajectories close 
to each other during their evolution. 



In analogy with standard definition (A2), we introduce 
a critical parameter 



A= lim (nt)- 1 Vln^, 



k=l 



(35) 



with do = |u>o(f)|- 
Now: 



(i) two particles start with identical initial conditions; the 
first one evolves according to geodesic equations (1)- 
(3), while the second one according to approximative 
equations (Eqs. (20)-(22), for example); 

(ii) i is determined as a time interval during which the two 
trajectories remain close to each other and the initial 
separation do is an arbitrary small number (typically, 
t is of the order of the Keplerian orbital period at the 
initial radius); 

(iii) Equation ( |35|) is evaluated and it is determined, 
numerically, whether convergence has been reached 
for large n with a pre-determined accuracy. The fi- 
nal value, Af, characterizes the rate of divergence of 
the two trajectories with the same initial conditions: 
Af > corresponds to separation increasing exponen- 
tially while Af — * — oo corresponds to a polynomial 
(i.e. much slower) increase of the separation. 

In order to estimate the accuracy of the approximative 
equations, one needs to follow a large number of trajecto- 
ries with randomly chosen initial conditions. By averaging 
over Af corresponding to different initial conditions, one 
obtains a value which shows whether the approximative 
equations describe, on the whole, the motion of test mat- 
ter with good precision. Regions where the mean value is 
positive, (Af) > 0, indicate that the approximation is not 
acceptable. The process is described in more detail below 
where we calculate (A f ) for Eqs. (11)-(13) and (20)-(22) 
as examples. 

4-2. An example: Test of the pseudo-Kerr equations 

We submitted our approximative equations to the above- 
described test on (Af). Relative position of the two test 
particles is determined by difference equations: 



Si 1 



J Korr 



approx' 



''Kerr 



approx' 



(36) 



with x l = {r(i), 0(f), 0(f)}. The suffix "Kerr" indicates 
that the trajectory is a geodesic in the Kerr space- 
time while the suffix "approx" corresponds to pseudo- 
Newtonian equations. We should stress in this place that 
comparing processes in different spacetimes is a serious 
problem in general relativity. Here, however, we do not 
compare two relativistic spacetimes and the situation is 
quite different: though a proper physical justification gives 
an additional appeal to any pseudo-Newtonian model and 
we have therefore emphasized also a physical content of 
our approach (in Sect. |^), in the final stage one mainly 
asks whether the test particles in the model field in some 
coordinates move along trajectories that are sufficiently 
close to certain geodesies of a given relativistic field in 
some particular coordinates. It is only necessary to de- 
fine the way of correspondence of the initial conditions. In 
our test, we consider as counterparts the particles which 
start from a given position (r,0) with given (specific) con- 
stants of motion (e = E/m,l — $/m). Of course, it is 
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possible that we would obtain a better fit with some other 
choice, e.g. if the "corresponding" particles had the same 
initial velocities with respect to some (corresponding) lo- 
cal frames. 




1.5 - 



1 10 100 1000 

Number of rescalings 



Fig. 1. This graph illustrates how the critical parameter A os- 
cillates at first and then settles to a positive value after k « f 3 
rescalingSjindicating the chaotic-type increase of separation dt 
(cf. Eq. (p5|)). Two different examples of trajectories are shown. 



(a) (b) (b) 




I / i 



Fig. 2. fsocurves of the mean terminal value of the critical pa- 
rameter, (Af), as a function of the particles' specific energy 
e and their specific axial angular-momentum I. Three pan- 
els correspond to (a) Eqs. (20)-(22) with a = M, (b) Eqs. 
(20)-(22) with a = 0, and, for comparison, (c) Eqs. (11)-(13) 
with a = M (scalar, Keres-Israel model). Positive (Af) indi- 
cates chaotic-type behaviour (see the text for details). 



We have integrated difference Eqs. ([56]) using the 
Bulirsch-Stoer scheme (Press et al. 1994). Initial condi- 
tions cover the parameter space of 

r > r+ ee M + y/M 2 -a 2 , < < §tt, e < e max . (37) 

Typically, we set e max ~ 1.1. As expected, most of the 
particles with e max 3> 1 escape quickly to large distances 
(r ^> r+) where both the exact and approximative equa- 
tions give identical results. It is thus relevant to investigate 
trajectories with lower energies which often make a num- 
ber of revolutions around the black hole. Some of these 
trajectories tend to diverge in a chaotic-type manner, as 
we will see below. 

Each run, characterized by the value of a/M, resulted 
in about 10 7 values of Af. For each set of initial condi- 
tions (^) within the run, the separation of the two cor- 
responding trajectories was periodically rescaled down to 
a pre-determined value after a fixed interval i. We have 
typically chosen i = 2M which is comparable to the Kep- 
lerian orbital period near the horizon. The integration was 
terminated when one of the following conditions had been 
satisfied: (i) A converged to a finite positive value Af (nu- 
merically, we checked that the relative change of A during 
the last 30 rescalings did not exceed 0.5%); (ii) the up- 
per limit of 10 4 on the number of rescalings was exceeded, 
while A was oscillating or decreasing monotonically to neg- 
ative values; (iii) the trajectory was captured by a black 
hole, r < r+. The case (iii) was excluded from further 
consideration. 

Figure [j] illustrates the evolution of A for two sets of 
initial conditions which both fall under item (i) above. Ap- 
parently, both cases shown there converge to positive Af, 
although the numbers of rescalings were different. In this 
respect, it is interesting to note that the separation of 
nearby geodesies in the Kerr spacetime never increases 
exponentially. This fact is a consequence of the existence 
of the fourth (Carter's) constant of motion, as discussed 
in Karas & Vokrouhlicky (1992). 

Next we illustrate the mean value (Af ) as a function of 
the constants of motion, e and I. Graphs corresponding to 
Eqs. (20)— (22) are shown for an extremely rotating black 
hole (Fig. |^a) and a non-rotating black hole (Fig. ||b). For 
comparison, the Keres-Israel model (11)— (13) with a = M 
was also included (Fig. ||c). One can locate easily the re- 
gions of positive (Af). In general these regions are bound to 
small e and I which correspond to trajectories that plunge 
close to the horizon. In the case of a non-rotating hole the 
graph is symmetrical about I = because both the Schw- 
arzschild spacetime and the adopted pseudo-Newtonian 
field with a — are spherically symmetric. Distortion in- 
troduced by rotation is visible in the graphs (a) and (c). 
Comparing these two graphs one can verify that model 
(19) is superior to ( jTfj| ) (cf. the positive values of Af in 
Fig. ||c, indicating that most of the corresponding trajec- 
tories separate rapidly from each other). 
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5. Conclusions 



We have described a method to test the overall accuracy 
of the pseudo-Kerr models of the relativistic gravitational 
field and we applied our approach to two particular po- n , 

tentials, given by Eqs. ( [To| ) and (19). One can conclude A = lim (nfp 1 ^ In — 
that the corresponding equations of motion, (11)— (13) and fe=l 

(20)-(22) respectively, do not provide an acceptable ap- 
proximation to exact Kerr geodesic equations in the part 
of the (e, Z)-plane where (A f ) > 0, and vice-versa. This 
statement does not apply strictly to all trajectories (with 
different initial conditions) because our graphs are given in References 
terms of mean values. We suggest that a check of this type 
should always be employed when some particular form of 
approximative equations is formulated. For example, cal- 
culation of the spectra of accretion discs requires integra- 
tion of a large number of photon trajectories; the overall 
accuracy in the description of the gravitational field is 
then important. 



rescaling). One can show (Benettin 1984) that the Lya- 
punov characteristic exponent corresponding to the origi- 
nal wo is given by 



(A2) 



independently of the value of t. In addition, A = A 
almost all Wq. Again, one can conveniently set d 



max 
1. 



for 
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A. Appendix 

The Lyapunov characteristic exponent A of two nearby 
trajectories, £ and £ + w, is defined by 



Chen X.-M., Igumen- 



1996, ApJ 461, 565 



165 



.4 



X(£,w) = lim i _1 ln 

t— too 



d(t) 



d(0) 



(Al) 



where £(t) is a solution of the equations of motion in the 
phase space, w(t) is a connecting vector, d{t) = \w\ is 
its length in the phase space and d(0) is an initial sep- 
aration (the value of d(0) must be small; otherwise it is 
arbitrary and one usually scales d(0) to unity). The Kerr 
spacetime being stationary, we can measure the separa- 
tions on the t = const surfaces. One defines A max as the 
maximum value of A with respect to variations of w and 
characterizes the chaoticity of the system by the value of 
Amax- Positive values of A max indicate that neighbouring 
trajectories diverge exponentially in the course of time 
while negative A max corresponds only to polynomial di- 
vergence. The above definitions have been originally in- 
troduced within the framework of non-relativistic systems 
but they are directly applicable also to stationary systems 
in general relativity. 

The maximum Lyapunov characteristic exponent is 
frequently determined numerically. This approach requires 
a careful choice of the integration scheme. In order to keep 
computational errors under control, one lets the trajecto- 
ries evolve for a short interval of time, i, after which w 
is rescaled back to unity (one denotes dk = \wk-i(t)\ the 
norm of the connecting vector at the moment of the k-th 
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